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Abstract 

A new cell- vertex scheme is outlined for solving for the flow about a delta wing with Moo > 1. Embedded regions 
of mesh-refinement allow solutions to be obtained which have much higher resolution than those achieved to date. 
Effects of mesh-refinement and artificial viscosity on the solutions are studied, to determine at what point leading-edge 
vortex solutions are “grid-converged.” A macroscale and a microscale for the size of the vortex are defined, and it 
is shown that the macroscale (which includes the wing surface properties) is converged on a moderately refined grid, 
while the microscale is very sensitive to grid spacing. The level of numerical diffusion in the core of the vortex is found 
to be substantial. Comparisons with experiment are made for two cases which have transonic cross-flow velocities. 


Introduction 

Compressible flow past sharp-edged delta wings can lead to a variety of flow patterns. Stanbrook and Squire [1] 
originally postulated the classification of flow patterns as a function of angle of attack normal to the leading-edge 

c*jv = tan” 1 (tan a sec A) 

and Mach number normal to the leading-edge 

Mjv = Mqo \/l — sin 2 A cos 2 a . 

Their data suggested a boundary near M/v = 1. To the left of the boundary ( M jy < 1) the leading-edge is swept 
behind the Mach cone and the flow is separated at the leading-edge. To the right (Mjv > 1) the flow is attached, 
with a Prandtl-Meyer expansion at the leading-edge. Szodruch and Peake [2] and Miller and Wood [3] carried the 
classification further, defining several regimes in the (Mjv,»at) plane. Miller and Wood [3], for example, produced the 
chart of flow regimes diagrammed in Figure 1. Vorropoulos and Wendt [4] postulated the possibility of leading-edge 
vortex flow with a reverse cross-flow shock under the vortex (See Figure 2), based on evidence from their laser- 
doppler results. This flow regime was at low normal Mach number and angle of attack. It is exactly this richness of 
flow phenomena that makes the compressible leading-edge vortex problem such an exciting one. The possibility of 
finding new complex flow topologies for this geometrically simple problem is a challenge, both for experimentalists 
and theoreticians. A good deal of numerical research has been carried out for this class of problems, as outlined below. 

Solutions of the Euler equations using centered-diflerence conservative schemes have been carried out by Rizzi [5], 
Rizzi and Eriksson [6], Murman et al [7] and Newsome [8], among others. Upwind conservative schemes have been used 
by Newsome and Thomas [9] and Chakravarthy and Ota [10]. Characteristic-based non-conservative schemes have 
been used by Marconi [11]. All of these give reasonable estimates for the pressure coefficient on the wing, although 
they tend to overpredict the suction peak due to the vortex. 

Some of these methods use the full three-dimensional form of the equations; others use a conical self-similarity 
assumption to reduce the problem to one with two independent spatial variables. The latter approach has justification 
through the experimental evidence described above: classification of the flows in the two-parameter (M/v, oat) space 
suggests conical self- similarity. Three-dimensional solutions for these problems suffer from a lack of resolution (Rizzi 
[5] has reported that at least 10 6 grid points are necessary to resolve the flow details), and even the conical solutions 
are not grid-converged [12]. To help understand and solve this problem, a solution method must be developed that 
allows a higher level of resolution than has been achieved previously. 
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The solution method described in this paper is a new cell- vertex scheme that allows for multiple levels of embedded 
grids. Cell-vertex schemes for logically rectangular meshes have been developed by Ni [13], Hall [14] and Siclari [15]. 
The embedded-grid scheme is adopted to provide good resolution in the leading-edge vortex region of the flow without 
the necessity of many grid points in the rest of the flow. The cell- vertex approach is adopted for its second-order 
accuracy, even at embedding interfaces and regions of highly stretched grid. Recent studies by Shapiro and Murman 
[16] have shown that a cell-vertex finite-volume scheme is identical to a finite elment method with constant test 
functions, and that its accuracy is equal to a Galerkin method for all practical purposes. The scheme as a whole 
allows highly accurate solutions to be obtained efficiently. Results from the scheme are used to show effects of grid 
refinement on leading-edge vortex solutions, and to determine the level of refinement necessary to resolve flow details. 


Governing Equations 


The three-dimensional Euler equations may be written in vector form as 
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where h 0 is the stagnation enthalpy, 

The ideal gas law 

p = (f-l)p^E- 

closes the set of equations. Introducing the conical variables 
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and assuming conical self-similarity (that the solution is independent of r), the Euler equations become 
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where 


u = u *f Y)V + C w 
v = v — rju 
w = w — (u 

and 

K = x/l + »? 2 +C 2 - 

These conical Euler equations may be expressed in terms of the state vector U, the Cartesian flux vector F and the 
conical flux vectors G and H as 

r <9U 8G dH 



The unsteady terms have been included so that an iterative procedure may be employed to reach a conically self-similar 
steady-state. The equations are solved on the unit sphere by setting r = 1. 

The physical boundary conditions consist of a no-flux condition at the wing, free-stream conditions outside the 
bow shock, and a Kutta condition at the leading-edge of the wing. Implementation of these boundary conditions will 
be discussed in the next section. 


Solution Procedure 

The solution scheme used is a finite- volume, multi-stage scheme in which the state variables are stored at the 
nodes. It allows for regions of local mesh refinement. The grid generation, spatial discretization, added artificial 
viscosity, temporal discretization and boundary and interface procedures are described below, along with the data 
structure that underlies the algorithm. 


Grid Generation 

For the cases presented here, the grid generation is carried out by a Joukowsky transformation. The ( 77 , C) plane 
is mapped to a complex x plane, in which the wing becomes a circle, by the transformation 

2 (§-A) 

»? + »< = X + tan 2 L 

2 * 

where A is the leading-edge sweep of the wing. In the \ plane, i ^constant lines are equiangularly spaced rays 
emanating from the origin and j ^constant lines are concentric rings. Grid points are generated along rays in the \ 
plane by an exponential stretching. This procedure yields near-conformal grids with good resolution near the wing. 
The embedded regions are generated by sub-dividing each cell in the region. This is done by adding a node on each 
face of the cell to be divided, and one in the center of the cell. Bilinear interpolations are used to find the coordinates 
of the new points. A typical grid is shown in Figure 3. Indexing the grid in an (ij) manner is no longer convenient; 
the data structure chosen for the scheme is described below. Although in the work presented here the embedded 
regions are determined a priori, the formulation and code are directly extendable to adaptative determination of the 
embedded regions. 


Discretization of Equations 

The finite volume discretization of the partial differential equations is formulated by integrating the conical Euler 
equations over a cell. This gives 
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Using Gauss’ theorem and the mean value theorem, this may be rewritten as 


r d U f r - * 

A — — + <f> Gd<-H<fy 
K dt Jqsi L 
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where an overbar denotes a cell-average and A is the cell area. 

The line integral of the fluxes is carried out by a trapezoidal integration about the cell, i.e. 


/ [cdC — Hd»j = [g (^1 + ^2) (»?2 ~ ^i) ~ g + ^2) (C 2 — Ci) 

faces *- 


where the subscripts denote the two nodes that define the face, ordered so that the integral is carried out in a 
counter-clockwise sense. \ 

The cell-average of the source term is calculat ed by a veraging the source term at the four nodes defining the cell. 
This gives the residual at the center of the cell, dU/dt. This cell-average of the residual is distributed to the four 
nodes defining the cell with a simple (4 weighting. 
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The artificial viscosity is a blend of a nonlinear second-difference and a linear fourth- difference. It is constructed 
from a weighted Laplacian, L y and an unweighted Laplacian squared, L 2 , respectively. The weighting function W is 
a normalized Laplacian of the pressure, 


>V = 


L(p)/p 

\\L(P)/P\\oo 


where || ||oo denotes the Loo norm, so that 0 < W < 1. This weighting is chosen for the efficient capturing of shocks 
and vortices. It causes the second-difference term to be first-order in high-gradient regions and small elsewhere; the 
fourth-difference term is third-order everywhere. Thus the artificial viscosity is given by 


D(IJ) = e 2 D 2 (U) - 6 4 D 4 (U) 

= e 2 I(W,U)-€ 4 £ 2 (U) 

where U is a modified state vector, with the energy term pE replaced by ph 0) so that the discrete equations permit a 
solution with constant total enthalpy. The Lapiacians are calculated in a way that automatically enforces conservation 
at physical and numerical boundaries. 

The temporal discretization is a four-stage scheme with coefficients ( 4 , 3 ,j,l)- The time step is calculated on a 
local basis to satisfy the CFL criterion for stability. 


Boundary Conditions 

There are boundary conditions to be enforced at the physical boundaries 

1. No flux through the wing 

2. Free-stream conditions upstream of the bow-shock 

3. Kutta condition at the leading-edges 

and at the numerical boundaries 

1. no flow through the symmetry plane 

2. consistency and conservation at the embedding interfaces 

The wing boundary condition is met by retaining only the pressure terms in the flux calculation on cell faces which 
abut the wing. No pressure extrapolation is necessary since the state variables are stored at the nodes. The outer 
boundary condition is implemented by ensuring that the outer boundary of the domain is outside the bow shock 
and enforcing free-stream conditions there. The Kutta condition is enforced implicitly by ensuring that the artificial 
viscosity is present near the leading-edges. This is done by arbitrarily setting the pressure switch W to one at several 
nodes in the vicinity of the leading-edges. Numerical experience demonstrates that this is adequate to ensure smooth 
separation [12]. 

The symmetry plane is introduced so that the flow past a wing at zero yaw may be solved on a half-plane. The 
condition to be enforced at the symmetry line is that the through-flow velocity, v ) is zero. This is implemented by 
setting v to zero initially and zeroing the distributed flux for the y-direction momentum equation at each iteration. 
Interfaces between two levels of grid refinement lead to cells of more than four nodes and faces. The embedding 
philosophy used here restricts the grid to a series of nested rectangles in which there is refinement by a factor of two 
in each direction at the interfaces. With this restriction, only one type of interface cell is generated — a cell with five 
faces and five nodes. These cells are the coarse cells which abut the embedding interface. This makes the interface 
treatment very simple. The flux integration for these cells is carried out exactly as above, with five faces for the 
cell instead of the usual four. The distribution step is modified, however, so that the “extra” node is updated. The 
distribution used is diagrammed in Figure 4. 


Data Structure 

Because of the embedded regions of the grid, it is inconvenient to index the grid in a typical (ij) fashion. Instead, 
nodes, cells, and faces are singly indexed, and pointer vectors are constructed with the necessary information for 
interrelating them. The pointer vectors used for the scheme are: 
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1. A cell-to-node pointer, relating a cell to the four corner nodes that define it 

2. A cell-to-face pointer, relating a cell to its four faces 

3. A face-to-node pointer, relating a face to the two nodes defining it 

4. A node-to-node pointer, relating a node to its (at most) four neighbors 

Each pointer has a specific use. The cell-to-node pointer is used in the flux distribution and the damping calculation. 
In each of these steps, cells are visited, and changes in the cell are distributed to the nodes of the cell. The cell-to-face 
pointer is used in the flux summation step — each cell is visited, and the fluxes from the four faces defining the cell 
are summed. The face-to-node pointer is used in the calculation of the face fluxes — each face is visited, and the 
coordinates and the values of the flux vectors at the nodes defining the face are used to carry out the trapezoidal 
integration. The node-to-node pointer is used to calculate grid metrics and cell areas, and to extrapolate quantities 
(such as the time step and the weighting function W) to the boundaries and interfaces. 

In addition to these, pointers which contain boundary information are necessary. They are 

L A vector containing the index of each face on the wing 

2. A vector containing the index of each node on the wing 

3. A vector containing the index of each node on the symmetry boundary 

4. A vector containing the index of each node on the far-field boundary 

5. A vector containing the index of each node at which the weighting factor for the damping is to be overridden 

6. A vector containing the index of each coarse cell which abuts an embedding interface 

7. A vector relating each of these interface cells to its “extra” node 

8. A vector relating each of these interface cells to its “extra” face 

9. A vector containing the orientation of each of these “extra” faces (i.e. which side of the cell the face is on) 

The boundary conditions at physical and numerical boundaries are enforced using these pointers. 


Effects of Numerical Parameters on Solutions 

The embedded mesh procedure outlined above forms a practical method to carry out a careful study of the effects 
of mesh refinement on leading-edge vortex solutions to the Euler equations. This study is carried out for a 75° swept 
wing at 10° angle of attack in a Mach 1.1 stream. The cross-flow streamlines for this case are shown in Figure 5. The 
vortex is the dominant feature of the flow. The cross-flow stagnation points at the leeward and windward symmetry 
points and on the outboard portion of the windward side of the wing are also evident. The cross-flow is transonic for 
this case. 

The grid is varied from an equivalent global resolution of 64 x 64 to 256 x 256. The results for the coarsest 
grid are shown in Figures 6-9. Figure 6 shows the grid, with refinement on the leeward side on the wing. Figure 7 
shows contours of the pressure coefficient. The compression on the leeward side of the wing, the expansion about the 
leading-edge and the low-pressure region in the core of the vortex may be seen. Figure 8 presents contours of the 
cross-flow Mach number. The large gradients of cross-flow Mach number in the core and feeding-sheet may be seen. 
Figure 9 presents contours of total pressure loss (1 — po/Poo)- There is a large loss in total pressure in the vortex, 
reaching a value of 34% in the core. 

Figures 10-13 are the results on a grid with equivalent 128 x 128 resolution. The gradients of pressure in the 
core (Figure 11) are better resolved, and a cross-flow shock is apparent underneath the vortex. The cross-flow Mach 
number contours are shown in Figure 12. Again the gradients are much better resolved, and the shock may clearly 
be seen. The contours of total pressure loss are shown in Figure 13. The lossy region appears smaller on this grid, 
but the level of loss has not gone down. Here it is 38%. 

A further grid refinement is shown in Figures 14-17. The pressure contours are shown in Figure 15. The shock 
is slightly sharper than on the previous grid and the suction in the core slightly better resolved, but otherwise the 
solution has not changed appreciably. The cross-flow Mach number contours (Figure 16) show the same trends. The 
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total pressure loss distribution has changed fairly drastically (Figure 17). The roll-up of the sheet into a core of 
distributed total pressure loss may be seen. The level of the loss is relatively unchanged, however, with a core loss of 
35%. 

The loss levels and two scales for the sizes of the vortices in the different cases are tabulated in Table 1. “Core 
sizes” are determined by summing the area of all cells in which the total pressure loss (1 — Po/POco) is greater than 10%. 
They are normalized by the core size for the 128 x 128 case. “Vortex width” and “Vortex height” are determined by 
measuring the width and height of the 2% contour at the center of the vortex. They are also normalized by the values 
for the 128 x 128 case. The “core size”, a microscale for the vortex, changes dramatically with the grid refinement. 
The “vortex width” and “vortex height” are nearly converged on the 128 x 128 grid, however. 

A similar table may be compiled for the effects of artificial viscosity on the loss. The values presented below are 
for a 128 x 128 grid and are normalized by the case for which c 2 = 0.003 and e 4 = 0.0010. Again the microscale for the 
vortex changes with the level of artificial viscosity, while the core loss level and the macroscale for the vortex remain 
approximately constant. The level of artificial viscosity does not greatly affect flow variables outside the core. 

It is interesting to note that, although the magnitude of the total pressure loss is independent of the level of 
artificial viscosity, it does depend on the form of the artificial viscosity. This dependence on the form of the viscosity 
was discovered in an attempt to lower the damping on the continuity equation. Figure 18 shows the total pressure 
loss contours for the M 0 0 = 1.1, ot = 10°, A = 75° case on an equivalent 128 x 128 grid where the second- difference 
damping terms in the continuity equation have been removed. For this case, the loss in the core is 27%. 

The effect that zeroing the continuity damping has on the form of the damping for the system of equations may 
be seen most easily by analyzing a one-dimensional steady model problem. The equivalent continuity and momentum 
equations are, including the second-difference damping, 
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Comparison with the one-dimensional steady Navier-Stokes equation shows that this equation has two spurious terms: 
the term u(e pu —e p )d 2 p/dx 2 } due to the difference in the two damping coefficients, and the cross-term 2e pu dp/dx du/dx. 
The first term is a spurious diffusive term; the second is a spurious convective term. 

If e pu = e p , the first spurious term drops out. This suggests a way to zero the continuity damping while retaining 
the proper form of damping for the system of equations. By subtracting these terms from the momentum and energy 
equations, the system will have the same form as it would in the case c pu = e p . Figure 19 shows the total pressure 
loss contours with the damping altered in this way. The core total pressure loss is 38%, as it was in the case without 
lowering the continuity equation damping. 

An estimate of the level of artificial viscosity in the vortex core may be determined by calculating an “equivalent 
Reynolds number” for the computation. In regions in which the second-difference viscosity dominates the fourth- 
difference viscosity, the damping terms for the conical Euler equations have the form 


D (u) ~ € 2 Z(w,u) 


~ e 2 WL ( 0 ) 

~ 4e 2 WA 2 V 2 U 


where A is an estimate of the grid spacing, and the factor of four arises from the non-normalized stencil used for the 
Laplacian operator L. Since this damping operator is multiplied by the CFL number in the temporal updating, the 
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actual damping terms added to the conical Euler equations have the form 

D (tr) ~ 4e 2 WAV 2 U 

where A is taken to be the average of the lengths of the faces emanating from the node. Thus an equivalent Reynolds 
number can be defined as 

- — = 4e 2 WA 

FL(Z e q 

Contours of the log of the equivalent Reynolds number are shown in Figures 20-22 for the three different levels of grid 
refinement. The equivalent Reynolds number increases with grid refinement, reaching values of 10 3 in the core for the 
256 x 256 case. The equivalent Reynolds number is surprisingly high at the leading-edge, due to the small size of the 
cells, but surprisingly low in the vortex core. Care must clearly be taken in computing solutions to the Navier-Stokes 
equations for leading-edge vortices in high Reynolds number flow. 

This study has shown that the level of total pressure loss in computed vortex cores depends on the form of 
the diffusive terms, but not on the magnitude of the diffusive terms. Analytical evidence to this effect has been 
demonstrated by the authors [17]. Numerous calculations performed by the authors show that the magnitude of total 
pressure loss depends strongly on the aerodynamic and geometric parameters [12]. Surface pressure coefficients and 
other variables outside of the core should converge with the macroscale for the vortex. For high-Reynolds number 
attributes of the vortex on a scale smaller than the macroscale, however, the grid resolution must be such that the 
equivalent Reynolds number due to truncation error and artificial viscosity is substantially higher than the Reynolds 
number of interest. 


Comparison with Experiment 

In this section, two cases are compared with experimental data. Both are 75° swept wings. They are 

1. Moo = 1-95, a = 25°, /? = 0° 

2. Moo = 1.70, a= 12°, /?= 8° 

Computed results for the first case are shown in Figures 23 and 24. Figure 23 shows the cross-flow streamlines for 
this case and Figure 24 shows the total pressure loss contours. Three cross-flow shocks are apparent in the vicinity 
of the vortex, and the vortex takes on an almost triangular cross-section. Figures 25 and 26 show a comparison of 
computed and measured pitot pressures. The pitot pressures compare extremely well in the vortex and above it; 
below it, viscous effects, which are not modeled in the code, dominate. Figure 27 shows a comparison of computed 
and measured flow angles along the leeward symmetry line. The agreement is surprisingly good, particularly near the 
wing. 

The second case is shown in Figures 28-33. This case was run on a grid with 256 x 128 equivalent refinement, 
shown in Figure 28. The cross-flow streamline plot (Figure 29) shows two vortices: a long, narrow vortex on the port 
side and a nearly axisymmetric vortex on the starboard side of the wing. There is a saddle point on the windward 
side of the wing near the port leading-edge, and one on the leeward side near the center of the wing. All the cross-flow 
streamlines to port of the two saddle points converge in the port vortex; all others converge in the starboard vortex. 
The total pressure loss contours (Figure 30) show the difference in strength and shape of the two vortices. The port 
vortex is stronger, with a core total pressure loss of 66%. The starboard vortex has a core loss of 43%. There is also 
evidence of a weak shock above the port vortex. The pressure coefficient contours (Figure 31) also show evidence 
of this shock. The pressure on the wing (Figure 32) compares well with experiment. Both vortices are predicted 
too far outboard, and the suction peaks are overpredicted. The solution does capture the unusual double expansion- 
recompression pattern under the port vortex, however. This unusual pattern is more clearly seen in the cross-flow 
Mach number distribution (Figure 33). Under the port vortex, the flow expands to a cross-flow Mach number of 1.6, 
recompresses to 1.2, then expands again, reaching a value of 1.3, then recompresses outboard to the port leading-edge. 
The cross-flow Mach number reaches a value of 1.7 above the port vortex, and goes slightly supersonic above and 
below the starboard vortex. 
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Conclusions 


A new cell-vertex scheme which allows for nested regions of embedded grid has been presented. The scheme has 
been used to solve several cases of compressible, supersonic flow past a delta wing. The grid-embedding capacity 
of the code has been used to study the effects of resolution on leading-edge vortex solutions. A macroscale and a 
microscale for the vortex have been proposed, and it has been shown that the macroscale converges on a moderately 
refined grid, while the microscale does not. It has been shown that even on a relatively fine grid, the level of artificial 
viscosity in the core of the vortex is substantial. Thus, if details inside the core are desired, it is necessary to use 
a very fine grid in that region to ensure that the equivalent Reynolds number due to truncation error and artificial 
viscosity is substantially higher than the Reynolds number of interest. 

The scheme has been used to compute two flows for which experimental data existed. Measurements of pitot 
pressure from the first experiment were modeled well, as was the flow angularity. Measured surface pressure coefficients 
from the second experiment were modeled fairly well, with discrepancies due to viscous effects. 
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Equivalent Grid 

Loss Level 

Core Size 

Vortex Width 

Vortex Height 

64 x 64 

33.9% 

2.55 

1.39 

1.43 

128 x 128 

37.6% 

1.00 

1.00 

1.00 

256 x 256 

35.1% 

0.26 

0.96 

1.05 


Table 1: Effects of grid refinement on loss 


*2 

^4 

Loss Level 

Core Size 

Vortex Width 

Vortex Height 

0.003 

0.0010 

37.4% 

1.00 

1.00 

1.00 

0.003 

0.0003 

37.6% 

0.32 

1.10 

1.13 

0.003 

0.0030 

38.2% 

1.74 

1.08 

1.02 

0.001 

0.0010 

37.2% 

0.71 

0.97 

1.01 

0.010 

0.0010 

37.6% 

1.67 

1.11 

1.08 


Table 2: Effects of artificial viscosity level on loss 
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Figure 3: Embedded Grid — Equivalent 128 x 128 Refinement 
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Figure 4: Interface Distribution Scheme 
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Afoo = 110 a = 10° 0 = 0 ° A =75° 128x128 
Cross-flow streamlines 
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Figure 5: Cross-Flow Streamlines 
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Figure 6: Grid — 64 x 64 Resolution 
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Moo = 1.10 a = 10° p = 0° A = 75° 64 x 64 
Pressure coefficient 
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Figure 7: Pressure Coefficient — 64 x 64 Resolution 
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Figure 8: Cross-flow Mach Number — 64 x 64 Resolution 
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Figure 9: Total Pressure Loss — 64 x 64 Resolution 
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Figure 10: Grid — 128 x 128 Resolution 
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Figure 11: Pressure Coefficient — 128 x 128 Resolution 
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Figure 12: Cross-flow Mach Number — 128 x 128 Resolution 
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Figure 13: Total Pressure Loss — 128 x 128 Resolution 
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Figure 14: Grid — 256 x 256 Resolution 
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Figure 15: Pressure Coefficient — 256 x 256 Resolution 
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Figure 16: Cross-flow Mach Number — 256 x 256 Resolution 
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Figure 17: Total Pressure Loss — 256 x 256 Resolution 
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Figure 18: Total Pressure Loss — Zero second-difference continuity damping 
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Figure 19: Total Pressure Loss — Zero second-difference continuity damping 
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Figure 20: Equivalent Reynolds Number — 64 x 64 refinement 
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Figure 21: Equivalent Reynolds Number — 128 x 128 refinement 
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Figure 22: Equivalent Reynolds Number — 256 x 256 refinement 
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Figure 27: Flow Angularity Comparison — M oo = 1.95, a = 25°, A = 75° 
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Figure 28: Grid — Moo = 1.7, a = 12°, A — 75°, /? — 8° 
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